We looked for combination therapies of neglected and misused antibiotics. Rifampin and polymyxin B combination stood out as a promising approach against scpectrum of clinical isolates. This repository contains all the data and code to reproduce the analysis. For details and full captions of figures, consult the publication.
We studied the potential of the rifampicin-polymyxin B combination against intra- and extracellular forms of bacteria: three P. aeruginosa strains, two clinical isolates of A. baumannii, E. coli, and K. pneumoniae.
dat = loadDdiData('input/dat/rds/ddi_all_strains.rds')
# Response surface analysis
# rs = getRS(dat, 'rs')
# rs = readRDS('input/dat/tmp/rs.rds')
#
# open3d()
# lapply(names(rs), sav3D)
#getDdiCounts(readRDS('input/dat/tmp/rs.rds')) %>%
# group_by(Strain) %>%
# ggplot(aes(Cond, Syn)) +
# geom_col(fill = 'grey70', col = 1) +
# facet_wrap(~Name, nrow = 3, dir = 'v') +
# labs(y = '# of synergistic concentrations', x = '') +
# theme(strip.text.x = element_text(size = 12))
#ggsave(filename='output/fig/SynergyCountAcid.svg', bg='white', h=7.5, w=7.5)
theme_set(theme_pdf())
sub = subset(dat, strain == 'ATCC27853' & (d1 == 0 | d2 == 0))
p1 = function() pltDR(effect ~ c1, subset(sub, c2==0), xlab = 'RIF, mg/L', zero=30)
p2 = function() pltDR(effect ~ c2, subset(sub, c1==0), col='orange', xlab = 'PMB, mg/L')
# pH
foo = fread("input/dat/raw/growth_and_pH.tsv")
foo = foo[time_h >= 0 & overnight_medium == 'pH7.4', .(
od = gmu(od620),
od.ci = gci(od620),
ph = mean(pH),
ph.ci = ci(pH)
), .(growth_medium, time_h)]
limit_of_detection = 0.03
foo[, od := replace(od, od < limit_of_detection, limit_of_detection)]
p0 = foo %>%
mutate(
growth_medium = relevel(factor(growth_medium), 'pH7.4'),
time_h = ifelse(time_h > 20, 20, time_h),
NULL
) %>%
ggplot(aes(time_h, od, shape = growth_medium)) +
scale_x_continuous(breaks = c(0, 5, 10, 20), labels = c(seq(0, 10, 5), 24), limits = c(0, 20)) +
scale_shape_manual(name = '', values = c(19, 22, 15), labels = c('pH 7.4', 'pH 5.5',
'pH 5.5 buffered')) +
theme(legend.position = 'top')
p3 = p0 +
scale_y_continuous(trans='log2', breaks = c(0.01, 0.03, 0.1, 0.3, 1, 3,
10)) +
geom_point(size=.point_size, col='grey20') +
geom_pointrange(aes(ymin = od*od.ci, ymax = od/od.ci)) +
geom_line(col='grey20') +
coord_cartesian(ylim = c(0.01, 10)) +
labs(x = 'Time, h', y = expression('OD'[620]))
p4 = p0 %+% aes(y = ph) +
scale_y_continuous(position = 'right') +
geom_point(size=.point_size, col=cbPalette[3]) +
geom_line(col=cbPalette[3]) +
geom_pointrange(aes(ymin = ph-ph.ci, ymax = ph + ph.ci), col=cbPalette[3]) +
theme(legend.position = 'none',
axis.title.y.right=element_text(color = cbPalette[3]),
axis.text.y.right=element_text(color = cbPalette[3]),
axis.ticks.y.right=element_line(color = cbPalette[3])
) +
labs(x = '', y = 'pH')
foo = fread('input/dat/raw/charcoal.csv')
labs = c('RIF, mg/L', 'PMB, mg/L',
rep(expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')'), 2)
)
tmpFun = function (i, df) {
par(mar=c(3.1, 6.1, 0, 0.7), mgp = c(1.7, 0.5, 0), tcl=-0.2)
if(i == 1) bar = subset(df, ab == 'rif')
else bar = subset(df, ab == 'pmb' | (ab == 'rif' & conc == 0))
# fit char vs no-char but share Emin and EC50 (argument pmodels)
fit = drm(response~conc, data=bar, CURVE,
pmodels=data.frame(CURVE,1,1),
fct=LL2.4( fixed=c(1, NA, NA, NA)))
plot(fit, type="average", cex.axis=1.25, cex.lab=1.25, lwd=2,
col='grey30', pch = c(1, 19), lty = 1, asp=1,
xlab = labs[i], ylab = labs[i+2], broken=T,
ylim=c(-4,4), legendText=c("control", "charcoal")
)
plot(fit, type="confidence", add=T, legend=F,
col=c('grey70', 'grey10')
)
}
p5 = function() tmpFun(1, foo)
p6 = function() tmpFun(2, foo)
aligned = align_plots(p3, p4, align = 'hv', axis = 'tblr')
plot_grid(
ggdraw(p1), ggdraw(p2), ggdraw(aligned[[1]]) + draw_plot(aligned[[2]]),
plot_grid(ggdraw(p5), ggdraw(p6), nrow=2),
labels = 'AUTO', align = 'v'
)
# ggsave('output/fig/SFig_mth_dr_pH_charcoal.svg', width=10, height=9)
# schematic for PD pars
tmp = data.table(
#Conc = c(0, 0.1, 0.3, 1, 3, 10, 30, 100),
Conc = c(0, 0.01, 0.03, 0.1, 0.3, 1, 3, 10),
Response = c(3, 3, 2.8, 1.5, -0.5, -2.5, -2.8, -3.2)
)
tmp[, Response := Response + 1]
# svg('output/fig/SFig_PDpar.svg', width = 3.5, height = 3.5)
#oldpar = par()
par(tcl = -0.3, mar=c(4.1, 4.1, 1.5, 0.5), mgp=c(1.5, 0.75, 0))
plot(
(mod = drm(Response ~ Conc, data = tmp,
fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
lowerl = c(-2.5, NA, NA),
robust = 'median'
)),
broken = T, pch = 19,
col = 'grey40',
ylim = c(-2.5, 4.7), xlim = c(0, 50),
xlab = c('Conc, mg/L'),
ylab = expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')'),
asp = 1
)
x1 = ED(mod, 0, type = 'absolute', logBase = exp(1))[1]
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:0 0.4249298 0.0092374
x2 = ED(mod, 50, logBase = exp(1))[1]
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 0.233915 0.005085
y2 = predict(mod, data.frame(Conc = x2))
segments(x1, 0, x1, -100, lty = 2)
segments(0.0001, 0, x1, 0, lty = 2)
segments(x2, y2, x2, -100, lty = 2)
segments(0.0001, y2, x2, y2, lty = 2)
abline(h = coef(mod)[2], lty = 2)
abline(h = coef(mod)[1], lty = 2)
text(x1+0.3, y = 0 + 0.1, expression(C[s]))
text(x2+0.3, y = y2 + 0.2, expression(EC[50]))
text(x = 0.002, coef(mod)[2] + 0.3, expression(E[min]))
text(x = 0.002, coef(mod)[1] + 0.4, expression(E[max]))
# dev.off()
# add Tn screen's PA14 (Liberati), and pmbR strains PA947 and PA1292
foo = rbind(
# pa14
fread('input/dat/raw/pa14_mut_pmb_rif_cfu.csv')[mut == 0] %>%
.[, `:=` (genus = 'Pseudomonas',
strain = 'PA14 Liberati',
mut = NULL)],
# pmbR
fread('input/dat/raw/pmb_resistant_isolates.csv') %>%
.[, genus := 'Pseudomonas']
) %>%
.[, `:=` (
c1 = d1*mic1,
c2 = d2*mic2)]
# transform
dat_fit = rbind(foo, dat, fill = T)[c1 == 0 | c2 == 0,
][, dose := ifelse(d1 == 0, c2, c1)
][, ab := ifelse(d1 == 0, 2, 1)
][, .(data = list(.SD)), .(strain, cond)]
# fit four parameter logistic regression
dat_fit[, fit := map(data, ~drm(effect ~ dose, ab, data = .x,
fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
pmodels = data.frame(ab, 1, ab),
lowerl = c(-4.5, NA, NA)))]
# get parameters
dat_fit[, pars := map(fit, ~getPDPar(.x))]
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 1.0936 0.7459 1.6033
## e:1:0 3017.5829 1432.1220 6358.2618
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 199.80084 143.78167 277.64579
## e:2:0 1.09877 0.85522 1.41167
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 4195.4112 844.9468 20831.4600
## e:2:0 1.4129 1.1031 1.8096
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 136.099 94.119 196.802
## e:2:0 108.833 77.678 152.485
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 116.627 78.836 172.534
## e:2:0 100.888 71.106 143.146
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 2.19824 1.73299 2.78839
## e:2:0 1.05512 0.74755 1.48923
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 0.54794 0.41895 0.71664
## e:1:0 2.54461 2.01152 3.21898
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 1.23558 0.90005 1.69619
## e:1:0 12.27513 9.50191 15.85774
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 0.72807 0.50460 1.05050
## e:1:0 8.19634 6.45160 10.41292
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 8.8108 5.2898 14.6755
## e:1:0 29.9697 20.8806 43.0152
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 2.9223 1.9854 4.3015
## e:1:0 25.5587 19.0791 34.2388
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 1.08515 0.83477 1.41063
## e:1:0 17.68410 13.82126 22.62653
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 3.0869 1.4710 6.4777
## e:1:0 9.1434 7.0217 11.9061
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 2.7579 1.6050 4.7391
## e:1:0 6.2385 4.5268 8.5975
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 0.71158 0.43400 1.16668
## e:2:0 4.78692 3.03655 7.54625
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 0.61538 0.37753 1.00308
## e:2:0 11.63563 7.65409 17.68830
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 10.31467 8.03876 13.23493
## e:2:0 0.52524 0.40478 0.68156
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 21.90092 17.02990 28.16518
## e:2:0 0.58875 0.45974 0.75395
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 34.814 26.015 46.590
## e:2:0 212.935 151.138 300.000
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 55.495 42.649 72.210
## e:2:0 24.544 18.867 31.930
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 22.7675 16.2560 31.8872
## e:2:0 2.0841 1.4483 2.9992
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 18.17315 12.63799 26.13261
## e:2:0 0.36616 0.23469 0.57127
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 6.273396 3.311455 11.884652
## e:2:0 0.071390 0.039096 0.130359
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 3.89931 2.70443 5.62209
## e:2:0 0.50220 0.38022 0.66331
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 5.68227 3.61600 8.92927
## e:2:0 1.27999 0.91009 1.80022
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 15.6419 11.5468 21.1894
## e:2:0 1.1580 0.8129 1.6495
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 34.12972 26.34817 44.20944
## e:2:0 0.56378 0.45798 0.69403
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 63.2718 45.7219 87.5581
## e:2:0 6.4509 4.8977 8.4967
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 78.80048 58.25397 106.59387
## e:2:0 1.15853 0.89306 1.50291
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 19.0752 14.1575 25.7011
## e:2:0 3.9380 2.8444 5.4521
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:1:0 38.15049 27.24067 53.42967
## e:2:0 0.92200 0.66327 1.28166
##
## Estimated effective doses
##
## Estimate Lower Upper
## e:2:0 0.60558 0.38290 0.95774
## e:1:0 27.16329 14.77923 49.92439
dat_fit = dat_fit[
dat_fit[, data[[1]][1, .(genus, mic1, mic2)], .(strain, cond)],
on = c('strain', 'cond')]
# reformat
dat_fit = unnest(dat_fit, pars) %>% as.data.table
dat_fit = getPdReadyForPlotting(dat_fit, dat)
p1 = dat_fit[!grepl('C', par)] %>%
.[lwr < -10, c('lwr', 'upr', 'value') := NA] %>%
ggplot(aes(value, strain, color = genus)) +
geom_vline(xintercept = 0, col = 'white') +
geom_pointrange(aes(xmin = lwr, xmax = upr), size = 1) +
geom_text(aes(label = round(signif(value, 2), 1)), col = 1, size = 2.5) +
facet_grid(cond ~ par, labeller = label_parsed,
scales = 'free', space = 'free_y') +
geom_blank(aes(x = min)) +
geom_blank(aes(x = max)) +
scale_color_manual(values = cbPalette[4:1], guide = guide_legend(
label.theme = element_text(angle = 0, face = "italic"), reverse = T)) +
scale_fill_manual(values = cbPalette[4:1]) +
labs(y = '', color = '') +
theme_html() +
theme(legend.position = 'top')
p2 = p1 %+%
mutate(dat_fit[grepl('EC', par)],
value = ifelse(value > 512, 512, value),
mic = ifelse(mic > 512, 512, mic)) +
geom_point(aes(mic, col = genus), size = 4.5, pch = 0) +
geom_text(aes(x = mic, label = round(mic, 1)), size = 2.5) +
scale_x_continuous(trans = 'log10', limits = c(0.03, 512),
labels = function(x) ifelse(x == 0, "0", x)) +
facet_grid(cond ~ par, labeller = label_parsed,
scales = 'free', space = 'free_y') +
xlab('Conc, mg/L') +
theme(legend.position = 'none')
plot_grid(p1 + xlab(expression(log[10]*'('*over('CFU'['24h'], 'CFU'['0h '])*')')),
p2, rel_widths = c(3, 2), align = 'hv', axis = 'tb')
# ggsave('output/fig/SFig_mth_pd.svg', width=10, height=6)
dtk = fread('input/dat/raw/time-kill.csv') %>%
remNonCharcoal() %>%
addNewVariablesToTimeKill()
# pltTimeKill(dtk)
#ggsave('output/fig/SFig_TimeKill.svg', width=15, height=5)
# # Print out the concentrations used
# unique(dtk[time_h > 0], by = c('broth', 'group')) %>%
# .[, .(broth, group, d_1, d_2)] %>%
# .[order(group)]
To account for the synergy in molecular terms—beyond a nonspecific increase in membrane permeability by polymyxin B—we turned to chemical genetics (Brochado and Typas, 2013). Working with ordered PA14 transposon library (Liberati et al., 2006), we derived a growth measure for monotherapies and combinations using colony opacity (Kritikos et al., 2017). To account for plate-to-plate variation, the opacity was multiplicatively corrected. This results in zero-centering of the Bliss scores, which were derived next. The significance of difference from zero Bliss score, for any mutant, was estimated by a T-test (5 biological replicates) and corrected for multiple testing (Benjamini-Hochberg).
tmp = chem.gen.res %>%
mutate(t.test.q.value = ifelse(t.test.q.value < 0.001, 0.001,
t.test.q.value)) %>%
mutate(showname = ifelse(
(t.test.q.value < 0.05 & (gene.name != '' | gene.name.to.show %in% c('66480', '26590', '43270', '02150'))) | t.test.q.value < 0.01 ,
T, F))
# tmp[t.test.q.value < 0.05 | t.test.q.value.gene < 0.05, .N, media]
# # media N
# # 1: LB 65
# # 2: LBpH5.5 39
# # only 3 are shared between the media
# foo = tmp[media == 'LB' & (t.test.q.value < 0.05 | t.test.q.value.gene < 0.05), Tn.mutant.id]
# bar = tmp[media == 'LBpH5.5' & (t.test.q.value < 0.05 | t.test.q.value.gene < 0.05), Tn.mutant.id]
# tmp[Tn.mutant.id %in% intersect(foo, bar)]
getLocus = function (dat, condition) {
unique(dat[eval(parse(text=condition))], by = c('locus'))[, locus]
}
by_gene = getLocus(tmp, 't.test.q.value.gene < 0.05')
by_mut = getLocus(tmp, 't.test.q.value < 0.05')
by_both = getLocus(tmp, 't.test.q.value < 0.05 | t.test.q.value.gene < 0.05')
by_gene_only = setdiff(by_gene, by_mut)
# # 70 hits
# length(by_both)
# # 53 from mutants
# intersect(by_mut, by_both) %>% length()
# # 17 additional from genes
# setdiff(by_gene, by_mut) %>% length()
# # 31 shared between mutants and genes
# intersect(by_gene, by_mut) %>% length()
Of the 70 hits, 53 are from mutant level of analysis and 17 come from gene level of analysis. 31 are shared between mutants and genes i.e. hits regardless if mutant or gene level of analysis is used.
(p0 = filter(tmp, media == 'LB') %>%
{
ggplot(., aes(epsilon.median.mutant.condition, -log10(t.test.q.value), col = t.test.q.value < 0.05)) +
geom_point(size = 2.5) +
scale_color_manual(values = c('grey70', 'grey20')) +
ggrepel::geom_text_repel(data= . %>% .[(showname)],
aes(label=gene.name.to.show), size = 6,
max.overlaps = Inf, box.padding=0.75, seed=3) +
lims(x = c(-0.8, 0.6), y = c(0, 3)) +
labs(x = 'Shift in Bliss score', y = expression('log'[10]*'(p-value)')) +
theme(aspect.ratio = 1, legend.position = 'none')
})
# ggsave('output/fig/ChemGen_LB.svg', width=5, height=5)
p0 %+% filter(tmp, media == 'LBpH5.5')
# ggsave('output/fig/ChemGen_LBpH5.5.svg', width=5, height=5)
| LB | LB pH 5.5 |
|---|---|
getGenes = function (dat) {
out = list()
for(i in unique(dat$media))
out[[i]] = unique(dat[t.test.q.value < 0.05 & media == i], by='locus') %>%
.[, (gene.name.to.show)]
return(out)
}
foo = getGenes(chem.gen.res)
#vtable = gplots::venn(foo)
With the following, we bring some biological knowledge into the analysis. This will get us at the level of processes/compartments as opposed to individual genes. We will focus on LB beacuse we have more data from there which is also more reliable (75% of the unique hits come from LB; there is less variance). In addition, our results suggest, the effect of pH on synergy is weak in PA14.
We did Gene Set Enrichment Analysis (GSEA) using GO terms from pseudomonas.com website and Kologorov-Smirnov testing for statistical significance estimation. Although most common approach, it has been critizised for example here.
library('topGO')
tab = getTopGO(chem.gen.res)
tmpFun = function (dat, name='CC') {
dat[[name]][1:10, c('GO.ID', 'Term', 'raw.p.value')] %>%
mutate(Pvalue = round(as.numeric(raw.p.value), 4)) %>%
dplyr::select(-raw.p.value) %>%
mutate(`GO.ID` = sub('GO:', '', `GO.ID`))
}
tmpFun(tab) %>%
knitr::kable()
| GO.ID | Term | Pvalue |
|---|---|---|
| 0008076 | voltage-gated potassium channel complex | 0.0096 |
| 0016021 | integral component of membrane | 0.0132 |
| 0055052 | ATP-binding cassette (ABC) transporter complex, substrate-binding subunit-containing | 0.0193 |
| 0005694 | chromosome | 0.0362 |
| 0033573 | high-affinity iron permease complex | 0.0530 |
| 0009289 | pilus | 0.1160 |
| 0005839 | proteasome core complex | 0.1221 |
| 0005960 | glycine cleavage complex | 0.1396 |
| 0030257 | type III protein secretion system complex | 0.1701 |
| 0005615 | extracellular space | 0.1717 |
# knitr::kable(format='latex', booktabs=T) %>%
# kableExtra::save_kable('output/fig/table_cell_component.png', font='Verdana')
# biological processes
tmpFun(tab, 'BP') %>%
knitr::kable()
| GO.ID | Term | Pvalue |
|---|---|---|
| 0009236 | cobalamin biosynthetic process | 0.0086 |
| 0009116 | nucleoside metabolic process | 0.0090 |
| 1900727 | osmoregulated periplasmic glucan biosynthetic process | 0.0166 |
| 0019700 | organic phosphonate catabolic process | 0.0173 |
| 0070475 | rRNA base methylation | 0.0179 |
| 0044249 | cellular biosynthetic process | 0.0260 |
| 0019354 | siroheme biosynthetic process | 0.0298 |
| 0006437 | tyrosyl-tRNA aminoacylation | 0.0346 |
| 0044262 | cellular carbohydrate metabolic process | 0.0357 |
| 0019748 | secondary metabolic process | 0.0370 |
# knitr::kable(format='latex', booktabs=T) %>%
# kableExtra::save_kable('output/fig/table_biological_process.png', font='Verdana')
Protein-protein interaction (PPI) analysis using STRING database. There was no data on PA14, so we will use PAO1 data to build and analyse the network onto which we then map PA14 orthologs.
library(igraph)
library(stringr)
hits = addAnnotation(chem.gen.res[t.test.q.value < 0.05, .(colony)])
hit_list = unique(hits[PAO1.ortholog != '', PAO1.ortholog])
# bring in the interactome
string_11_pa = 'input/dat/raw/pao1_ppi_stringDB_287.protein.links.v11.0.txt' %>%
fread() %>%
.[, `:=` (
combined_score = combined_score/1000,
#strip / trim species ID
protein1 = gsub("287.DR97_", "PA", protein1, fixed=TRUE),
protein2 = gsub("287.DR97_", "PA", protein2, fixed=TRUE)
)]
# create an igraph object
pa_ppi = graph.data.frame(string_11_pa, directed = FALSE)
# summary(pa_ppi)
V(pa_ppi)$degree = degree(pa_ppi)
pa_ppi = simplify(pa_ppi, remove.multiple=T, remove.loops=T, edge.attr.comb='mean')
# extract network for all the hits in the provided list
hit_ppi = induced.subgraph(pa_ppi, which(V(pa_ppi)$name %in% hit_list))
# store degree on the network (main hit network)
V(hit_ppi)$degree = degree(hit_ppi)
# Replace the PAO1 genes with corresponding PA14 ones
lookup = unique(hits[, .(PAO1.ortholog, gene=gene.name.to.show)], by='PAO1.ortholog')
V(hit_ppi)$name = left_join(data.table(PAO1.ortholog=V(hit_ppi)$name), lookup)[, gene]
# Communities ------------------------------
# Community detection based on edge betweenness (Newman-Girvan)
# pdf('output/pdf/graph.pdf', height=10, width=20)
par(mfrow=c(1, 1), mar=c(1,3,1,1), cex=0.6)
set.seed(1)
ceb = cluster_edge_betweenness(hit_ppi)
V(hit_ppi)$community = ceb$membership
cols = ins(cbPalette, 'white', c(5, 8:9, 11:14))
#svg('output/fig/SFig_ppi_cluster.svg', width=20, height=20)
#dendPlot(ceb, mode="hclust", palette = myCols, cex=1.5, cex.axis=1.5)
plot(hit_ppi, vertex.color=cols[V(hit_ppi)$community],
vertex.label.cex=2, edge.width=4, vertex.label.family="Helvetica",
vertex.label.dist=0, vertex.size=16,
vertex.frame.color='black', vertex.label.color='black'
)
#dev.off()
The major graph communities, using (Newman-Girvan’s edge betweenness):
We validate the sensitivity of identified candidate mutants in low throughput and in liquid LB medium at pH 7.4. Instead of factorial (i.e. checkerboard), we use a fixed ratio design Tallarida et al 1997.
dat = fread('input/dat/raw/pa14_mut_pmb_rif.csv') %>%
addDrugRatio() %>%
addAUC() %>%
addFitness() %>%
remDissimilarExperimentalConditions() %>%
addPA14MutAnnotation()
dat[, dose:= ifelse(cond_f %in% c('PMB', 'combo'), d1, d2)]
dat_lst = split(dat, dat$mut)
fit = lapply(dat_lst, fitPA14MutDR, curveid = cond_f)
# Loewe response surface ------------------------------
# rs = getPa14RS(dat)
# We'll use a snapshot of results from 2021-01-02
rsl = readRDS('input/dat/rds/2020-01-02_rs_pa14.rds') %>% lapply(., function(x) x$rsl)
foo = lapply(rsl, function(x) summary(x$maxR)$totals) %>% rbindlist(idcol='mut')
# Mutant numbering is potential source of confusion; hence I assign here
# manually: wild-type is coded as mut 0, fitting of mut 28 failed
foo[, mut := c(0:27, 29:45)]
#arrange(foo, Syn) %>% print(20)
pltNSynergiesOfValidatedMutants()
# ggsave('output/fig/validated_mutants_synergy_count.svg', w=5, h=10)
#svg('output/fig/SFig_45PA14MutDoseResponses.svg', 12, 16)
#plt45PA14MutantDR() # dose-responses
#dev.off()
#svg('output/fig/SFig_45PA14MutCompLoeweNull.svg', 12, 16)
# We'll use a snapshot of results from 2021-01-02
#plt45PA14MutantComparisonToLoeweNull('2020-01-02_rs_pa14')
#dev.off()
getValidHitTable() %>% addEcOrthologs() %>% .[,-'locus_pa'] %>%
knitr::kable()
| PA gene | EC ortholog | Synergy | Description |
|---|---|---|---|
| 02150 | none | Serine phosphatase RsbU, regulator of sigma subunit | |
| relA | relA | none | (p)ppGpp synthetase |
| 26590 | none | GntR family transcriptional regulator | |
| 66480 | less | Predicted ATPase involved in chromosome partitioning | |
| 56840 | less | Predicted thiol oxidoreductase | |
| lldP | glcA | less | L-lactate permease |
| 60490 | less | cytochrome c | |
| 03760 | more | sodium:solute symporter | |
| 43270 | selU | more | tRNA 2-selenouridine synthase |
| 51310 | more | Predicted redox protein, regulator of disulfide bond formation | |
| maiA | more | maleylacetoacetate isomerase | |
| 43670 | more | sensor/response regulator hybrid | |
| pmbA | pmbA | more | PmbA protein |
| 62230 | more | Predicted kinase |
# knitr::kable(format='latex', booktabs=T) %>%
# kableExtra::save_kable('output/fig/table_valid_hit.png')
getGoTermTable() %>%
knitr::kable()
| PA gene | Location | Process | Function |
|---|---|---|---|
| 02150 | membrane | signal transduction | catalytic |
| 03760 | membrane | transmembrane transport | transmembrane transporter |
| 26590 | NA | regulation of transcription; biosynthesis | catalytic; transcription factor |
| 43270 | NA | tRNA seleno-modification | transferase for selenium-containing groups |
| 43670 | membrane | signal transduction; phosphorylation | phosphorelay sensor kinase |
| 56840 | NA | NA | electron transfer; heme binding |
| 60490 | NA | NA | electron transfer; heme binding |
| lldP | membrane | lactate transport | lactate transmembrane transport |
| maiA | cytoplasm | aromatic amino acid metabolism | catalytic; protein binding |
| pmbA | NA | peptidoglycan biosynthesis; proteolysis | metallopeptidase |
| relA | NA | ppGpp metabolic process | NA |
#knitr::kable(format='latex', booktabs=T) %>%
#kableExtra::save_kable('output/fig/table_valid_hit_go.png', font='Verdana')
#rs = fread('input/dat/raw/pa14_mut_pmb_rif_cfu.csv') %>%
# prepForRS() %>%
# setDT() %>%
# addMutGeneNames() %>%
# getRS(., 'rs_6mut')
rs = readRDS('input/dat/tmp/rs_6mut.rds')
rsl = lapply(rs, function(x) lapply(x, function(y) y$rsl))
# lapply(names(rs), sav3D, xcap=10)
memb = readRDS('input/dat/rds/permeability.rds')
pltMemb(group_by(memb, PmbConcMgL, Strain))
#ggsave('./output/fig/PA14MembranePermeability.svg', w=15, h=5)
memb = readRDS('input/dat/rds/permeability.rds')
om = addAovDunnetMembrane(memb)
im = addAovDunnetMembrane(memb, membrane = 'IM')
getDunnetResultsMemb(om, lev = 1)
## $PmbConcMgL
## [1] 0.1
##
## $RifConcMgL
## [1] 0
##
## $Dunnett
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`wild-type`
## diff lwr.ci upr.ci pval
## 66480-wild-type -2.106667 -6.2881962 2.0748629 0.4809
## 02150-wild-type -3.770000 -7.9515296 0.4115296 0.0819 .
## relA-wild-type -1.123333 -5.3048629 3.0581962 0.8942
## 26590-wild-type -1.423333 -5.6048629 2.7581962 0.7821
## 43270-wild-type 3.376667 -0.8048629 7.5581962 0.1304
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## $RobustDunnett
## contrast estimate SE df t.ratio p.value
## 66480 - (wild-type) -2.11 1.92 12 -1.098 0.6889
## 02150 - (wild-type) -3.77 1.91 12 -1.974 0.2417
## relA - (wild-type) -1.12 3.01 12 -0.373 0.9773
## 26590 - (wild-type) -1.42 1.92 12 -0.743 0.8698
## 43270 - (wild-type) 3.38 1.92 12 1.758 0.3293
##
## P value adjustment: dunnettx method for 5 tests
getDunnetResultsMemb(im, lev = 1)
## $PmbConcMgL
## [1] 2
##
## $RifConcMgL
## [1] 0
##
## $Dunnett
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`wild-type`
## diff lwr.ci upr.ci pval
## 66480-wild-type -7.536667 -16.605398 1.5320642 0.1162
## 02150-wild-type -9.753333 -18.822064 -0.6846025 0.0341 *
## relA-wild-type -25.990000 -35.058731 -16.9212692 2.2e-05 ***
## 26590-wild-type -29.013333 -38.082064 -19.9446025 5.3e-07 ***
## 43270-wild-type 10.973333 1.904602 20.0420642 0.0169 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## $RobustDunnett
## contrast estimate SE df t.ratio p.value
## 66480 - (wild-type) -7.54 2.79 12 -2.704 0.0735
## 02150 - (wild-type) -9.75 5.71 12 -1.708 0.3523
## relA - (wild-type) -25.99 3.07 12 -8.471 <.0001
## 26590 - (wild-type) -29.01 2.14 12 -13.588 <.0001
## 43270 - (wild-type) 10.97 2.72 12 4.041 0.0069
##
## P value adjustment: dunnettx method for 5 tests
dat = 'input/dat/raw/pmb_resistant_isolates.csv' %>%
fread() %>%
setnames(., 'strain', 'mut') %>%
prepForRS(swapnames=T) %>%
setnames(., 'mut', 'strain') %>%
setDT()
# rs = getRS(dat, 'rs_pmbR_intra')
# lapply(names(rs), sav3D)
dat = loadDdiData('input/dat/rds/ddi_all_strains.rds')
dat_fit = anti_join(dat, dat[c1==0 & c2!=0])[, .(data = list(.SD)), .(strain, cond)]
dat_fit[, fit :=
map(data, function(x) {try(
drm(effect ~ d1, d2, data = x,
fct=LL2.4(names=c("h", "m", "b", "e"), fixed=c(1, NA, NA, NA)),
lowerl = c(-4.5, NA, NA))
)})
]
## Error in optim(startVec, opfct, hessian = TRUE, method = "L-BFGS-B", lower = lowerLimits, :
## non-finite value supplied by optim
## Error in drmOpt(opfct, opdfct1, startVecSc, optMethod, constrained, warnVal, :
## Convergence failed
idx = sapply(dat_fit$fit, function(x) class(x) != 'try-error')
dat_fit[idx, ed90 := map(fit, function(x) list(try(ED(x, 90, interval='fls', display=F) %>%
round(3))))]
dat_fit[idx, par90 := lapply(unlist(dat_fit[['ed90']], recursive=F), function(x) try(as.data.table(x, keep.rownames='level')))]
par90 = lapply(1:nrow(dat_fit[idx, ]),
function(x) dat_fit[idx, ][['par90']][[x]]) %>%
rbindlist(idcol='id') %>%
merge(dat_fit[, .(strain, cond)][idx][, id := 1:26], .) %>%
lowerCaseNamesDT() %>%
.[, level := gsub('e:|:90', '', level)] %>%
.[, effect:= 90]
set.seed(1)
res = par90 %>%
group_by(id) %>%
mutate(rel_IC = estimate/head(estimate, 1),
rel_IC = ifelse(rel_IC > 2, 2, rel_IC)
) %>%
mutate(level = ifelse(level < 1 & level != 0, 0.1, level)) %>%
filter(as.numeric(level) <= 1)
# res %>% group_by(level) %>% summarise(ec90=median(estimate), quantile(estimate, 0.2), ci=gci(estimate), lwr=ec90/ci, upr=ec90*ci) %>% dplyr::select(-ci)
# level ec90 lwr upr
# 1 0 7.41 6.34 8.65
# 2 0.1 3.20 2.32 4.41
# 3 1 2.18 0.966 4.92
# res %>% filter(level==0.1) %>% pull(estimate) %>% quantile(0.2)
# 20%
# 0.939
ggplot(res, aes(level, estimate)) +
geom_jitter(width=0.2, pch=21, aes(fill=cond), col='grey70') +
geom_boxplot(alpha=0.5, outlier.shape=NA) +
scale_y_continuous(name=expression('EC'[90]^'RIF'*', mg/L'),
trans='log2', breaks=c(16, 8, 4, 2, 1, 0.5, 0.25)) +
scale_x_discrete(name='PMB, xMIC', labels=c('0', '< 1', '1')) +
coord_cartesian(ylim = c(0.2, 16)) +
labs(fill = '') +
theme(aspect.ratio=1)
#ggsave('output/fig/rif_ec90_upon_pmb.svg', width=5, height=4, bg='white')